Structure of the magma plumbing system beneath Semisopochnoi Island (Aleutian Arc) inferred from seismic tomography

Semisopochnoi Island is a remote and little-studied volcanic island in the western part of the Aleutian Arc. The existence of several active volcanic centers and a 5000–7000-year-old large caldera makes this island an important site for volcanic hazard assessment in the Northern Pacific. Based on local seismicity data recorded by six permanent seismic stations, we created the seismic tomography model, including the 3D distributions of Vp, Vs, and Vp/Vs ratios to a depth of 10 km. This model provides the first geophysical insight into the interior structure of Semisopochnoi Island and sheds light on the processes in the magma plumbing system beneath all volcanic centers on the island. At depths of 5–10 km, we observed a columnar-shaped high Vp/Vs-ratio anomaly below the caldera in the central part of the island, which likely represents the steady magma conduit. This conduit is headed by a prominent high Vp/Vs-ratio anomaly located 3–5 km directly below the caldera, which represents the magma reservoir feeding Cerberus and other Holocene-aged volcanic centers on Semisopochnoi Island.

The high level of volcanic and seismic activity of the Aleutian Arc (Fig. 1a) may possibly pose serious problems for the entire North Pacific region 1,2 . The volcanic islands of the Aleutian Arc are characterized by the presence of traces of violent explosive eruptions in the recent geological past 3 . Potential eruptions and earthquakes in this area can pose a threat to people living on the islands and dense aviation routes 4 . To respond to this threat in a timely manner, networks of telemetric seismic stations operate permanently in most key Alaskan and Aleutian locations and continuously transmit data in real time to the office of the Alaska Volcano Observatory (AVO) 5 . These stations have been in operation for decades, so they provide useful information on the internal processes and structures of active volcanoes. Data from these networks have been used to create seismic tomography models for several Aleutian and Alaskan volcanoes, which has helped better understand the evolution of the magma sources [6][7][8][9][10][11] .
Besides natural hazard issues, geophysical studies of structures beneath volcanic areas are important for the exploration of geothermal resources, which are some of the most prospective sources of renewable carbon-free energy. Although most of the Aleutian Islands are inhabited, geothermal power plants built on them might also be used to generate hydrogen, which is in increasing demand due to the ecological challenges facing the world.
Here, we consider the remote and little-studied Semisopochnoi Island located in the western part of the Aleutian Arc. Episodic volcanic activity involving multiple vents makes this island an important site for the volcanic hazard assessment of the Western Aleutians. Furthermore, the existence of a large 5000-7000-year-old caldera 12 indicates that the volcanoes on Semisopochnoi Island could potentially produce catastrophic explosive eruptions with regional and, possibly, global impacts. Due to the difficulty of access, a limited number of experimental studies have been performed there. The first geological map of Semisopochnoi Island was compiled by Robert Coats and his team based on two expeditions conducted in the 1940s and the 1950s 13 12 report, which summarizes the results of these studies, is the main source of information for the short geological overview presented in the next section.
Most importantly for us is the permanent network of telemetric seismic stations that the Coombs' team deployed during the 2005 campaign, which has operated continuously since then. Using the arrival times of P-and S-waves from the local seismicity recorded by this network, we created a seismic tomography model of the crust beneath Semisopochnoi Island. As far as we know, this is the first geophysical image of the crustal structure beneath Semisopochnoi Island. The importance of this study is that it sheds light on the details of the magma plumbing systems of the Semisopochnoi Island volcanoes. Moreover, a comparison of the results of this study with those of other Aleutian volcanoes will provide a better understanding of the common principles of the Aleutian Arc's functioning arc volcanism. Data. In this study, we used the local seismicity P-and S-wave arrival times from January 1, 2013, to December 31, 2017, provided by the AVO. These data were recorded by six permanent Semisopochnoi Island stations, installed in September 2005 22 . The catalog contained information on 3410 events with 11,711 P-and 9803 S-wave arrival times (average of 6.3 picks per event). We used several criteria to select the data for the tomography, including (1) the distance from an event to the nearest station should be less than 20 km, (2) the total number of P-and S-wave picks per event should be equal to or larger than 4, and (3) the time residual after locating the events in the starting 1D velocity model should be smaller than 0.5 s. After applying these criteria, the number of data was significantly reduced to 1784 events (6142 P-and 5606 S-wave picks; 6.6 picks per event). www.nature.com/scientificreports/ The locations of the events after running the full tomography procedure are shown in map view and two vertical sections in Fig. 2. Note that some of events in the vertical sections appear to be above the surface, because we projected the events onto the profile plane from a certain band, where higher topography features than those on the profile may exist. The low number of stations used in the current study compared to most other tomography studies may have resulted in the poor source locations and the low quality of the tomographic images. However, other studies with similar station configurations have provided fair-quality results, such as eight stations on Avacha volcano 25 , eight stations on Nevado del Huila 26 , and four stations on Udina 27 . All of these studies underwent subsequent extensive testing that proved the reliability of the models. Following the same strategy, we undertook a number of tests to demonstrate the resolution limitations of our case. www.nature.com/scientificreports/ The distributions of the P and S wave velocities, as well as the Vp/Vs ratio, were calculated using the LOTOS code for passive-source tomography 28 , which was previously implemented in a number of volcanoes with similar data configurations. Details on the algorithm are presented in "Methods" section.

Results
The results of the synthetic tests allowed us to assess the spatial resolution of the model and the accuracy of the source locations. This was particularly important given the relatively low number of stations. For the synthetic modeling, we reproduced all the conditions of the experimental data processing, including the same calculation steps and controlling parameters. Here, we present two groups of tests aimed at separately assessing the vertical and horizontal resolutions. Synthetic travel times using the same source-receiver pairs used for the experimental dataset were calculated in a predefined synthetic model using bending ray-tracing. These travel times were perturbed by random noise, with an average deviation of 0.05 and 0.1 s for the P-and S-wave data, respectively. This provided a variance reduction similar to the inversion of the experimental data. We then "forgot" all the information about the true velocity model, source coordinates, and origin times, and reproduced all the experimental data processing steps, including finding the initial locations of the sources in the starting 1D model.
The results of three horizontal checkerboard tests with different sized anomalies are shown in Fig. 3. In all cases, the anomalies were defined as unlimited columns, and their size and the spacing between them are indicated in the figure captions. The anomaly amplitude was ± 7%, and it had the opposite signs for dVp and dVs. The recovered values of the Vp/Vs ratio are shown in two horizontal sections at depths of 2 and 7 km. The 3-km and 2.5-km anomalies were robustly resolved at both depth levels, whereas the 2-km anomalies appeared to be slightly smeared. This test indicated the minimum size of the anomalies that could be resolved in different parts of the study area, based on the existing data.
In another test, three synthetic models were defined along three vertical sections, the same as those used to present the main results in Fig. 4. In this case, the models consisted of two layers of alternating anomalies with the sign change interval at a depth of 3-5 km. In the direction along the section, the size of the anomalies was 1.5 km, and the empty spacing between the anomalies was 2 km. The thickness of the anomalies in the direction across the section was 6 km. Similar to the previous case, the ± 7% amplitude anomalies of the P-and S-wave velocities had opposite signs, which enabled contrasted variations of the Vp/Vs ratio. All these models were robustly recovered in the central part of the model, including the 3-5 km deep transition. However, instead of www.nature.com/scientificreports/ narrow anomalies separated by large interval of zero values, we observe some smearing, when the transitions between anomalies appear to be smooth. Similar effect might exist in the case of the experimental data inversion. The synthetic tests provided the opportunity to estimate the accuracy of the source locations. During the recovery procedure, the true source locations were presumed to be unknown; therefore, the difference between the final and true locations provided an adequate image of the real uncertainty of the source locations during the experimental data inversion. An example of source mislocations for the vertical checkerboard along Sect. 3 (shown in the right-hand column of Fig. 4) is presented in Fig. 5. After locations in the starting 1D model, the mean error was 0.65 km; however, after five iterations of the simultaneous inversions of the source and velocity parameters, the error was reduced to 0.47 km. As expected, the uncertainty of the source depth locations was larger than that of the horizontal coordinates. Note that despite such significant errors and outliers, the velocity model was correctly recovered.
The results of the experimental data inversion are presented in Figs. 6 and 7 in four horizontal and three vertical sections with the Vp and Vs anomalies and the Vp/Vs ratios, which were calculated by the division of the resulting absolute P-and S-wave velocities. The stability of this method for calculating Vp/Vs was confirmed by the series of synthetic tests described above. The final locations of the sources of this velocity model are shown in map view and vertical sections in Fig. 2. www.nature.com/scientificreports/

Discussion
The derived dVp, dVs, and, especially, Vp/Vs-ratio anomalies provided important information about the magma plumbing processes beneath the Semisopochnoi Island volcanoes. Based on experimental studies 29 and several previous tomographic studies of other volcanoes [8][9][10][11]25,30,31 , some common relationships between seismic parameters and rock properties can be proposed. For example, P-wave velocity is mostly sensitive to rock composition, whereas S-wave velocity is mostly affected by the presence of a liquid phase (e.g., melt and volatiles). Therefore, semi-molten magma stores and conduits are usually associated with anomalies with higher dVp, lower dVs, and a very high Vp/Vs ratio. For the most active volcanoes worldwide, where P-and S-wave velocities are available, such relationships are very helpful for imaging magma reservoirs and conduits 8-11,25,30,31. However, this pattern might not be valid at shallow depths, as there are too many different processes affecting the mechanical properties of the rocks, such as the interlayering of soft sediments and rigid lava flows, the penetration of meteoric water, and hydrothermal alteration. In this case, a separate interpretation of Vp and Vs appears to be more informative than consideration of the Vp/Vs ratio. Indeed, as shown in Fig. 6 in the sections at 0 and 2 km deep, most of the Holocene craters are located inside or on the border of the low-velocity Vp and Vs anomalies, but they are almost not expressed by the Vp/Vs ratios. A prominent low Vp and Vs anomaly is located in the center of the caldera, directly below Windy Cone. This anomaly may represent the accumulation area of the soft volcaniclastic rocks that filled the caldera and/or a hydrothermally altered zone. The Cerberus Cones are located on the border between the low-velocity anomaly and a high-velocity pattern, representing rigid, highly consolidated rock. Such a mechanically differential contact between media represents a weak zone, through which magma ascent is most plausible. An association of active volcanic vents with borders between high-and low-velocity zones is observed in many volcanoes (e.g., La Palma volcano 32 and Tolbachik volcano 33 ).
In the deeper sections at 4 and 7 km (Fig. 6), the general configuration of the anomalies appears considerably different compared to the shallow sections. Here, the Vp and Vs anomalies generally look inversely correlated. The positive Vp anomaly inside the caldera coexists with the strong negative Vs anomaly, which causes a very high Vp/Vs ratio, exceeding 2. These patterns are typical attributes of magma storage areas observed below many active volcanoes 25,30,31 . In all the vertical sections shown in Fig. 7, this anomaly is located between 3 and 5 km   www.nature.com/scientificreports/ This anomaly matches well with the results of the mechanical inverse modeling, indicating that the source of the magma that created the ~ 25-cm inflation was located at a depth of 5-10 km below the central part of the caldera 24 , exactly where we located the high Vp/Vs-ratio anomaly. The vertical sections in Fig. 7 show how the Holocene-aged Semisopochnoi Island vents were fed. In "2" section of the Vp/Vs ratio, the main 3-5 km deep reservoir has several inclined upgoing branches. The left branch ascends to the Ringworm and Threequarter monogenic cones, and the right branch ascends toward the Windy Cone. However, the elevated Vp/Vs-ratio anomaly was less prominent than the periphery branches below the Cerberus stratovolcano, where most of the historical activity was observed. An explanation for this could be that most of the Cerberus eruptions were mainly basaltic. Low-viscosity basaltic magma can propagate through narrow fractures; it does not require large conduits. Therefore, the magma migration paths beneath Cerberus cannot be resolved in our model. Another explanation could be the different sensitivities of the P-and S-wave velocities at shallow depths. They are both negative below Cerberus, but the existence of other near-surface processes prevents the origin of high Vp/Vs patterns.
Another important feature is observed beneath the Sugarloaf volcanic complex, another historically active volcanic center. In all the horizontal sections below the Sugarloaf Cones, a local high Vp anomaly coexists with a strong low Vs anomaly and a high Vp/Vs ratio, a clear attribute of semi-molten magma. As shown in vertical "3" section in Fig. 7, the vertical high Vp/Vs-ratio anomaly beneath the Sugarloaf Cones appears to be isolated from the shallow magma source beneath Cerberus. It is possible that they have a common source at a greater depth that cannot be resolved with the existing data. www.nature.com/scientificreports/ The structure obtained for Semisopochnoi Island in this study is consistent with a number of other tomography models derived from other Aleutian and Alaskan volcanic islands. For example, Koulakov et al. 's 11 modeling of the structure beneath Akutan Volcano indicated a single deep conduit up to 5 km below the surface that branches at shallower depths to form a number of local conduits connected to volcanic vents and centers of geothermal activity. Koulakov et al. 's 9,31 tomography modeling beneath Mount Spurr indicated a columnar high Vp/Vs-ratio anomaly from 4 to 5 km below the surface, which, at shallower depths, is replaced by a zone of low Vp/Vs. The existence of well-established single conduits in these cases might be explained by the steady state of the Aleutian subduction, which has remained in the same location for hundreds thousand or even millions of years. This would determine a common scenario for the development of arc volcanism that started with the long-term accumulation of basaltic material from a single conduit over millions of years. This finally produced a large shield volcano, creating a circular island. As the crust thickened, the composition of the magma became more silicic. At some point, a strong explosive eruption formed a large caldera. As the upper crust became more heterogeneous, several magma pathways were created that connected the deep conduit to the surface. As a result, volcanic eruptions of variable compositions occurred through multiple vents distributed across a large area, both inside and outside the caldera.

Conclusion
We have presented the first seismic tomography model for Semisopochnoi Island. It provides insight into the structure beneath the island and reveals details of the magma plumbing system beneath the island's Holoceneaged volcanoes, some of which are currently active. Because we used data from only six seismic stations, we performed a number of synthetic tests to prove the sufficient resolution of the tomography model.
In the deep part of the model, a columnar anomaly of high Vp, low Vs, and high Vp/Vs ratio revealed the location of the magma conduit, which was consistent with the location of the magma source estimated from the modeling of ground deformation during the 2014 seismic swarm 24 . Above this conduit, we observed a very high Vp/Vs-ratio anomaly reaching depths of 3-5 km located beneath the presently active Cerberus volcano, which likely represents its magma reservoir. This anomaly appears to be connected to other Holocene-aged volcanic centers located inside and outside the caldera. However, another active volcanic complex, Sugarloaf, seems to be isolated from the Cerberus magma reservoir and located above another vertical conduit.
This tomography model has potential for the exploration of geothermal energy sources on Semisopochnoi Island. The upper limit of the magma reservoir is located ~ 3 km beneath the central part of the caldera, so a 2-2.5 km deep borehole would likely provide a sufficient temperature gradient to enable the profit-making exploration of geothermal energy, which can then be used for hydrogen production.
An important methodological conclusion of this study is that the low number of stations available to us was sufficient for the satisfactory quality of the calculated seismic model. This opens up the opportunity to obtain interesting results for other areas with seismic small networks that were previously considered too small for seismic tomography inversions.

Methods
To undertake iterative simultaneous inversions of the velocity model and the source locations, we used LOTOS, passive-source tomography code 28 , which was used for the cases cited above and many other studies. We used the standard version of the LOTOS code, which has been described in detail in previous articles; therefore, here, we will only present the major steps with indications of the key controlling parameters.
The calculations start with determining the first approximations for the source locations using the grid search method and a simplified algorithm to calculate travel times along straight lines. To search for the best source locations, we successively refined the grid (i.e., 4 × 4 × 4 km, 1 × 1 × 1 km, and 0.1 × 0.1 × 0.1 km). We then updated the source determinations using a 3D ray-tracing algorithm based on the bending method 34 .
The 3D distributions of the P-and S-wave velocity anomalies were parameterized by a set of nodes installed in the study area according to the distribution of the ray paths. We set a minimum grid spacing of 1 km and 0.5 km in the horizontal and vertical directions, respectively. In both cases, the grid spacing is much smaller than the minimum size of anomalies that can be recovered with the existing data. Therefore, the resolution is determined by damping and smoothing, and not by grid spacing. The smaller vertical spacing is defined to make the grid more flexible to describe velocity variations near the surface, where significant relief variations exist. To decrease the effect of the grid geometry to the results, we performed a series of inversions in grids with different basic azimuths (i.e., 0°, 22°, 45°, and 66°) and then averaged the results in a single 3D model with a regularly spaced grid.
For each grid, we performed a simultaneous inversion for the P-and S-wave velocity anomalies (i.e., dVp and dVs) and the source corrections (i.e., three parameters for the coordinates and one for the origin time). The derived matrix was inverted using the least squares method LSQR algorithm 35,36 . To control the stability of the model, we used amplitude damping (i.e., 0.2 and 0.4 for the P and S models, respectively) and flattening (i.e., 0.4 and 1.9 for the P and S models, respectively). The values of these parameters were determined based on the results of the synthetic modeling.
After the inversion step, we updated the 3D velocity model and used it as a new reference model for the next iteration, which started with the relocation of the sources. We then calculated the new matrix and performed the inversion again. In total, we performed five iterations, an optimal compromise between the quality of the solution and the calculation time. In the final step, the average residual deviations in the L1 norm were reduced from 0.077 to 0.070 s for the P-waves and from 0.071 to 0.062 s for the S-wave data. The starting 1D model was defined at several depth levels and linearly interpolated between levels. We used a constant Vp/Vs ratio of 1.85. The P-wave velocities at each level were determined after several runs of the full inversion procedure by adjusting the average velocity according to the results of the previous model. The initial